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ABSTRACT 

Future high redshift 21 -cm experiments will suffer from a high degree of contamination, due 
both to astrophysical foregrounds and to non-astrophysical and instrumental effects. In order 
to reliably extract the cosmological signal from the observed data, it is essential to understand 
very well all data components and their influence on the extracted signal. Here we present 
simulated astrophysical foregrounds datacubes and discuss their possible statistical effects on 
the data. The foreground maps are produced assuming 5° x 5° windows that match those 
expected to be observed by the LOFAR Epoch-of-Reionization (EoR) key science project. 
We show that with the expected LOFAR-EoR sky and receiver noise levels, which amount to 
f« 52mK at 150 MHz after 400 hours of total observing time, a simple polynomial fit allows 
a statistical reconstruction of the signal. We also show that the polynomial fitting will work 
for maps with realistic yet idealised instrument response, i.e., a response that includes only 
a uniform uv coverage as a function of frequency and ignores many other uncertainties. Po- 
larized galactic synchrotron maps that include internal polarization and a number of Faraday 
screens along the line of sight are also simulated. The importance of these stems from the fact 
that the LOFAR instrument, in common with all current interferometric EoR experiments has 
an instrumentally polarized response. 

Key words: cosmology: theory, diffuse radiation, observation, radio lines: general, instru- 
mentation: interferometers, radio continuum: general 



1 INTRODUCTION 

The Epoch of Reionization (hereafter, EoR), which marks the end 
of the Universe's 'Dark Ages', is one of the least explored epochs 
in cosmic evolution. Currently, there are two main observational 
constraints on the EoR. The first is the sudden ju mp in the Lyman- 
Q op tical depth in the Gunn-Peterson troughs jOunn & Petersonl 
Il965[) observed in the Sloan Dig i tal Sky Survey quasar spec- 
tra iBeckeretal.' '200 1'; 'Fan et al. 1200 ll : IPentericci et"ai] 120021 : 



I Whit e et al.q 2003; Fan et al. 2006), marking a lower limit to the 
redshift at which the Universe became completely ionized. The 
second constraint comes from the fifth year WMAP satellite data 
on the temperature and polariz ation anisotropies of the cosmic mi- 
crowave background (CMB) ( Spergel et al.l2007l : |Page et al.l2007h 
which gives an integral constraint on the Thomson optical depth 
experienced by the CMB photons since the EoR. However, both 
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of these observational methods provide limited information on the 
reionization process. 



The redshifted 21 -cm hyperfine transition line of 
neutral hydrogen is the most promising and immedi- 
ately accessible method for probing the intergalactk 
medium (IGM) during reionization (e.g. Field 195^ 
Scott & Rees 1990; Kumar, S ubramanian. & Padmanabha 
1995; .Madau. Meiksin. & Reea .19971) . Recent years have wit- 
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nessed a flurry of theoretical activities to pre dict reionization 
sourc es and their impa c t on t he IGM (e.g. 'Barkana & Loeb| 
I2OOII; iLoeb & Barkarial I2OOII; ICiardi, Ferrara & W hite 200," 
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Bromm & Larsonll20o4 illiev etd] 
Thomas & Zaroubil l2007h . Mea 



Ciard i. Stoehr. & White 200 
120071 ; IZaroubTet aL qOO* * 
surements of the 21 -cm signal can al so help to constrain the 
cosmological parameters independently jMcOuinn et"ai]|2006h . 
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Future telescopes like LOFAkQ, MW^E 21Cm4!I and SKJ^ 
are being designed to study the redshifted 21 -cm signal from the 
EoR. A successful detection of this signal will help us derive the 
nature of the first sources and their impact on the surrounding IGM. 

Unfortunately however, the cosmological EoR signal is con- 
taminated by a slew of astrophysical and non-astrophysical com- 
ponents. Typically, the contamination level is orders of magnitude 
larger than the cosmological 21 -cm signal. Thus, the primary chal- 
lenge of the EoR observations will be the accurate modelling of 
the various data components - foregrounds, instrumental response, 
ionospheric disturbances, to name a few - which is essential to de- 
velop a robust signal extraction scheme. 

For the foregrounds, there are currently no available data in 
the 115-180 MHz frequency range and 4 arcmin resolution at 
high Galactic latitude that would allow accurate modelling of the 
LOFAR-EoR foregrounds. Therefore, one has to rely on the avail- 
able relevant data and extrapolate, based on theoretical arguments, 
into the frequ ency range and resolution obse rved by LOFAR. How- 
ever, recentlv lXii. Bharadwai. & Chengahirl |2008) used 153 MHz 
obervations with Giant Meter Wave Radio Telescope to character- 
ize the statistical properties - visibility correlation function - of the 
foregrounds. This paper focuses on simulating the galactic and ex- 
tragalactic foregrounds that dominate the sky at frequencies of in- 
terest for the LOFAR-EoR experiment (115-215 MHz). The main 
foreground components are: Galactic synchrotron emission from 
diffuse and localised sources. Galactic thermal (free-free) emis- 
sion and integrated emission from extragalactic sources (like radio 
galaxies and clusters). The dominant component of the foregrounds 
is the Galactic synchrotron emission (~70 per cent). The extra- 
galactic emission contrib utes ~ 27 per ce nt and Galactic free-free 
emission ~1 per cent ( Shaver et al. 1999). Although the difference 
between the mean amplitude of the EoR signal and the foregrounds 
is expected to be 4-5 orders of magnitude, an interferometer like 
LOFAR measures only the fluctuations which in this case are ex- 
pected to be different by 'only' three orders of magnitude. 

Various authors have stu died the foreg r ounds in the con- 
text the EoR measurements. IShaveretal.1 ( Il999l) have stud- 
ied the diffuse synchrotron and f ree-free emissio n from our 
Galaxy and ex tragalactic sources; iDi Matteo et al.l ^20021) and 
lOi Matteo. Ciar di. & Miniati ( 20045 have considered emission 
from unreso lved extragalactic sourc es at low radio fr equencies; and 
lOh & Mackl ( [2003,) and Coorav & Furlanettd j2004l) studied the ef- 
fect of free-free emission from extragalactic haloes. Over the years, 
several methods have been explored to filter out the foregrounds. 
Most of the methods rely on the relative smoothne ss in the fre 



quency of the foregrounds, with respect to the signal |Shaver et"al 
19991; jpi Matteo et alj|2002l; IZaldarriaga. Furlanetto, & Hemguisl 



2OO4I; [Morales. Bowmari.& Hewittl |2006| ; IWangetal.1 l200a 

Gles er. Nusser. & Benson 200?'). 

[SantpSjJI^ooray, & Knox (2005) have studied the foregrounds 
for the EoR experiment and their influence on the measurement 
of the 21-cm signal. In their multi-frequency analysis of the power 
spectra, they considered four types of foregrounds: Galactic diffuse 
synchrotron emission; Galactic free-free emission; extragalactic 
free-free emission; and extragalactic point sources. They showed 
that foregrounds cleaning is aided by the large scale angular cor- 
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relation, especially of the extragalactic point sources, which facili- 
tates signal extraction to a level suitable for the EoR experiments. 

The current study is part of the general effort undertaken by 
the LOFAR-EoR key science project to produce simulated data 
cubes. The pipeline under construction will simulate the LOFAR- 
EoR data cube that includes the simulated cosmological 21-cm 
signal, the galactic and extragalactic foregrounds, ionospheric ef- 
fects, radio frequency interferences (RFIs) and the instrumental re- 
sponse. These datacubes will be used to design the observational 
strategy and test our signal-processing methods. Our main concern 
in this paper is the simulation of the galactic and extragalactic fore- 
grounds. 

Recently, a study by iGleser. Nusser. & Bensonl ( l2007h has 
been conducted along lines similar to parts of the current paper. The 
authors test a certain signal extraction algorithm on simulated fore- 
grounds maps in which they take most of the relevant foregrounds 
into account. However, there ar e many important differences be- 
tween the two papers. First, in the lGleser. Nusser. & Bensonl i2007h 
study the assumption for the noise level in the LOFAR-EoR project, 
as well as the other experiments, is at least an order of magnitude 
too low. They assume 1 and 5 mK noise models whereas in reality 
the noise for the LOFAR-EoR experiment is about 50 mK. They 
also present a simplified model of the Galactic foregrounds that 
does not take into account all the spatial and frequency correlations 
of the Galactic diffuse synchrotron emission and underestimates 
that of the Galactic free-free emission, both of which are very im- 
portant. In contrast to them, we also present polarized maps and 
introduce the LOFAR instrumental response and noise in a realistic 
manner. 

In the foregrounds simulations presented in this paper we 
choose a different approach from previous groups, since our main 
aim is to produce the simulations that will be part of the LOFAR- 
EoR data pipeline. In this context our main aim is to produce fore- 
grounds maps in the angular and frequency range of the LOFAR- 
EoR experiment, i.e. 3D datacubes, and then use those simula- 
tions for testing the accuracy of removal of the foregrounds. Sec- 
tion|5]outlines the importance of the polarized character of the fore- 
grounds and how to model the Stokes I, Q, and U polarization maps 
of the Galactic synchrotron emission. Section[6]presents simulated 
instrumental effects of the LOFAR telescope and their influence on 
the foregrounds maps, and Section|7]discusses a method to extract 
the EoR signal from the foregrounds. The paper concludes with a 
discussion and outlook (Section[8j. 



2 THE COSMOLOGICAL 21-CM SIGNAL 

In radio astronomy, where the Rayleigh- Jeans law is applicable, 
the radiation intensity, is expressed in terms of the brightness 
temperature Tt, such that: 



-kTb, 



(1) 



where u is the frequency, c is the speed of light and k is Boltz- 
mann's constant. The predicted differential brightness tempera- 
ture deviation of the cosmological 21-cm sign al from the cosmic 
microwave backgroun d radiation is given by (|Fii3[l95i,[l95i; 
ICiardi & Madauli2003l) : 

Sn = 26mKxHi(l +5) ( 1 
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Here Ts is the spin temperature, xm is the neutral hydrogen frac- 
tion, 5 is the matter density contrast, Q.m and f^b are the mass and 
baryon density in units of the critical density and h = _ffo/10oE|- 

In his seminal papers. Field ( 1958, 1959) used the quasi-static 
approximation to calculate the spin temperature, Ts , as a weighted 
average of the CM B, kinetic and colour temperature jWouthuvsenl 
ll952ljFieldll958h : 

rp _ TcMB + VkinTkin + VaTa 
1 + Vkin + Va 

where Tcmb is the CMB temperature and and j/^ are the ki- 
netic and Lyman-Q coupling terms, respectively. We have assumed 
that the color temperature, Ta, is equal to Tkin- The kinetic cou- 
pling term increases with the kinetic temperature, whereas the ya 
coupling term is due to t he Lyman-g pumping, known a lso as the 
Wouthuysen-Field effect dWouthuvsenll 1952l : lFieI(3l 19581) . The two 
coupling terms are dominant under different conditions and in prin- 
ciple could be used to distinguish between ionization sources, e.g., 
between first stars, for which Lyman-Q pumping is dominant, vs. 
first mini-quasars f or which X-ray photons and therefore heating is 
dominant (see e .g., 'Nussei'^2005':'Kuhlen, Madau, & Montgomervl 
l200dlZaroubr et al. 2007; Thomas & Zaroubi 2007). 

The brightness temperature of the cosmological signal used 
in this study is produced from a dark-matter-only N-body simula- 
tion. This simulation is used to produce a cube of the cosmolog- 
ical signal, i.e., the density as a function of right asces ion, decli- 
nation, and redshift (for more details see lThomas et all in prepa- 
ration). Although Ts is calculated according to Eq.[3] we assume 
that Ts ^ TcM B ■ The reason for this assumption is that towards 
the redshifts of interest for the experiment (z =6-12), the abun- 
dance of Lya photons in the Universe is sufficient to co uple Ts to 
Tk which is obviously much greater than Tcmb ( Ciardi & Madaul 
[2003) . Hence from Eq.[2l follows the cosmological density and 
xhi- We further assume that along each sight-line the neutral frac- 
tion follows the function 1/ (1 + exp(z — Zreion)), where Zreion 

for each pixel (or line of sight) is set to 8.5±Sz=io and where 5^=10 
is the density contrast at redshift 10. We used this approach to ran- 
domize the reionization histories along different lines of sight while 
preserving the spatial correlations of the cosmological signals. In 
principle, this randomization could be drawn out of a Gaussian dis- 
tribution function. Redshift 10 here is an arbitrary choice. Zreion 
along each line of sight varies in accordance with the cosmological 
density along that line-of-sight at z=10 and has a variance of unity 
centred at 8.5. Fig.[T]shows the signal data cube that we use in order 
to test our foregrounds filtering procedure. 

Currently, a number of experiments (e.g., LOFAR, 2ICMA, 
MWA and SKA) are being designed to directly measure STt of the 
HI 2 1 -cm hyperfine line and probe the physics of the reionization 
process by observing the neutral fraction of the IGM as a function 
of redshift. In this study, we focus on predictions for LOFAR, but 
our conclusions could be easily applied to the other telescopes. 

TheLOFAR-EoFQkey project plans to measure the brightness 
fluctuations in the frequency range of 115-190 MHz, correspond- 
ing to redshift range 6-11.5 with spectral resolution of ~ 1 MHz 
and angular resolution of about ~ 4 arcmin. A more detailed dis- 
cription of the LOFAR array will be given later in the paper when 
the instrumental effects are discussed (Section|6f. 

5 We assume a ACDM Universe with = 0.04, = 0.26, = 
0.738 and i^o = 70.8 k ms-^Mpc"! 

^ For more information, see the LOFAR web site: www.Iofar.org and the 
LOFAR-EoR web site: www.astro.rug.nl/ LofarEoR. 
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Figure 1. Simulated EoR signal assuming an exponential form for the reion- 
ization history as decribed in the text and Ts 2> Tqmb- The simulation 
box is 100 Mpc (comoving) a side. The upper panel shows the differ- 
ential brightness temperature in a slice along the redshift/frequency direc- 
tion and another spatial direction. The lower panel shows the brightness 
temperature as a function of redshift/frequency along a certain sight line 
(the dashed line in the upper panel). The resolution along the frequency 
direction is 10 kHz. 

3 GALACTIC FOREGROUNDS 

The Galactic foregrounds have three main contributions. The first 
and largest component is the Galactic diffuse synchrotron emis- 
sion (GDSE), which is the dominant foreground component in the 
frequency range of the LOFAR-EoR experiment. The second com- 
ponent is radio synchrotron emission from discrete sources, mostly 
supernova remnants (SNRs). The third and last component is the 
free-free radio emission from diffuse ionized gas. This component 
is the weakest of the three, yet it still dominates over the cosmolog- 
ical component. Moreover, it has a different spectral dependence, 
making it very imporant in testing the signal extraction schemes 
that we have. In this section we describe how we simulate the con- 
tribution of each of these components to the total intensity. The 
polarized intensity simulations are described later on. 

3.1 Galactic diffuse synclirotron emission (GDSE) 

The GDSE originates from the interaction between the free elec- 
trons in the interstellar medium and the Galactic magnetic field. 
Therefore the observed GDSE intensity as a function of frequency, 
depends on the number density of emitting electrons, Ne, and 
the Galactic magnetic field component perpendicular to the line of 
sight, B±_ : 

I{v) ~ N,b'^2^'^'\~^'<-^^'^ (4) 

whe re 7 is the electron spectral energy distribution power law in- 
dex l |Pacholczvklll970l) . The intensity of the synchrotron emission 
as expressed in terms of the brightness temperature varies with po- 
sition and frequency and its spectrum is close to a featureless power 
law Tb ~ , where /3 is the brightness temperature spectral index, 
related to 7 by /3 = -(2 -|- (7 - 1) /2). 

Observationa l data that are relevant to th e LOFAR-EoR 
project are scarce. iLandecker & Wielebinskil ( ll970h have produced 
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an all sky map of the total intensity of the GSDE at low radio 
frequencies at 150 MHz with 5° resolution. The other Galactic 
survey releva nt to the L OFAR-EoR experiment is the 408 MHz 
survey of Hasla m et alj (|l982) with a resolution of 0.85° and of 
Reich & Reich (198^ at 1420 MHz with 0.95° resolution. In the 
Reich & Reich 419881) paper the authors also assume a smooth 



power law change in the intensity as a function of frequency which 
they calculate from their 1420 MHz and 408 MHz maps. 

At high Galactic latitudes the minimum brightness tempera- 
ture of the GDSE is about 20 K at 325 MHz with variations of 
the order of 2 p er cent on scales from 5-30 arcmin across the sky 
jde Bruvn et. a l 1998). At the same Galactic latitudes, the temper- 
ature spectral index [3 of the GDSE is abou t —2.55 at 100 MHz 
and steepens toward s higher frequencies (e.g. lReich & ReichI 19881 : 
IPlatania et al.l Il998h . Furthermore, the spectral index gradually 
changes with position on the sky. This change appears to be caused 
by a variation in the spectral index along the line of sight. An ap- 
propriate standard deviation in the power law index, a/j, in the fre- 
quency range 100- 200 MHz appears to be of the order of ~ 0.1 
( IShaveretalJl999l) . Recent data, collected around a galaxy cluster 
Abell 2255 using the WSRT telescope at 350 MHz, indicate that 
the rms of the brightness temperature at 3 arcmin resolution could 
be as low as 0.1-0.3 K (Pizzo and de Bruyn, private communica- 
tion). If extrapolated to 150 MHz this result implies that the rms in 
that region could be 1-2 K, which is an order of magnitude smaller 
than the low resolution data suggest. 

For the purpose of this paper we assume that the GDSE as a 
function of frequency is well approximated by a power law within 
the limited frequency range of 115-180 MHz. This is a central as- 
sumption in our simulation which is consistent with the general 
trend shown by the available data, namely that the change in the 
frequency power law index is gradual. The values we choose for 
the po wer law index are ba sed o n the high Galac t ic lati tude regions 
in the iHaslam eTZI ( Il982l) and iReich & Reiclj h988l) maps. The 
second assumption we make is that both the intensity and power 
law index of the GDSE can be spatially modelled as Gaussian ran- 
dom fields (GRFs). For the power spectrum of GRFs we assume a 
power law with 2D index ct — —2.7. The standard deviation of the 
GRFs is normalized to 0.4, assuming an angular scale correspond- 
ing roughly to the field of view (5°). This is consistent with the 
value adopted bvlTeg mark et alJ ( l2000l) . lGiardino et alj ^20021) and 
ISantos, Coorav, & Rnox t2005i) for the angular power spectrum in- 
dex a, where Ci ~ a varies from —2.4 to —3, and / is the 
harmonic number. 

In contrast to the previous authors I Tegmark et alj l2000l : 
iGiardino et"Z]|2002l : ISantos. Coorav. & Knoxll2005h whodirectlv 
used the angular and frequency power spectrum of the GDSE for 
their analysis, we simulate GDSE in four dimensions (three spatial 
and one frequency), produce maps at each frequency and then do 
our analysis on them. The four dimensional realisation approach 
has the added benefit of enabling us to account for the ampli- 
tude and temperature spectral index variations of the GDSE along 
the line of sight (^-coordinate). We obtain the final map of the 
GDSE at each frequency, i/, by integrating the GDSE amplitude 
(A{x, y, z, v)) along the z-coordinate: 



Th{x,y,u) = C / A{x,y,z,v)Az 



(5) 



where Tb{x, y, v) is the brightness temperature of the GDSE as a 
function of position and frequency and C is a normalization con- 
stant. A{x, y, z, v) is dimensionless and at each frequency is de- 
fined by power law: 
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Figure 2. Simulated map of total intensity emission of Galactic diffuse syn- 
chrotron emission with angular spectral index a = —2.7 and frequency 
spectral index /3 = —2.55. The angular size of the map is 5° X 5°, with 
~ 0.6' resolution. The colour bar represents the brightness temperature Ti, 
of the Galactic diffuse synchrotron emission in Kelvin at 120 MHz. 



A{x, y, z, v) = A{x, y, z, i^o) — 



!3{x ,y .z ,1^) 



(6) 



where i/q is the reference frequency at which the normalisation is 
done and f3{x, y, z, v) is the temperature spectral index as a func- 
tion of 3D position and frequency v. The power law index /3 has a 
weak frequency dependence, also as a power law. 

A{x, y, z,uo) and P{x, y, z, vo) of the GDSE at the normal- 
ization frequency uq are modelled spatially as two Gaussian ran- 
dom fields with 3D power law spectrum P{k) ~ fc*. Note that the 
absolute value of the 3D power law index 5 is |(5| = \a\ + 1 where 
a is the 2D power law index mentioned above. A{x,y, z, ^o) and 
P{x, y, z,uo) are normalized according to observations (the Galac- 
tic surveys mentioned above). 

For clarity, the steps we followed to produce the GDSE maps 
are listed below: 

(i) Generate the same 3D Gaussian random field for both A and 
f3. The assumption here is that both fields have a correlated spatial 
distribution, which is supported by visual inspection of the high 
Galactic latitude portions of the iReich & Reicj ( 11988 ^ maps. We 
have also explored the possibility that A and /? are independent; this 
has led to results very similar to the correlated case, and therefore 
we show only maps in which A and (3 are correlated. 

(ii) Normalize the mean and standard deviation of A{x, y, z, i/q) 
by integrating along the z direction and setting the mean and stan- 
dard deviation of Tt{x,y ,1^0) to match the observations (the Galac- 
tic surveys mentioned above). In other words we set the integra- 
tion constant C in Eq. |5] in a way that the properties of the field 
A{x,y,z,h'o) after integration match the observed properties of 
Tb{x,y,iyo). 
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Figure 3. Brightness temperature of Galactic diffuse synchrotron emission 
as a function of frequency, for five different lines of sight. Each line of sight 
has a slightly different power law index along the frequency direction as 
a result of the spatial and frequency variations in the temperature spectral 
index. 

(iii) Normalize the mean and standard deviation of 
(i{x, y, z, uq) according to observations. 

(iv) Use Eq.[6]to calculate A at each frequency. 

(v) Integrate along the ^-coordinate to get the two-dimensional 
maps of the GDSE brightness temperature at each frequency f 
(Eq.|5}. 

For the purpose of this paper we simulate the GDSE on 512^ 
grid, where the x,y plane corresponds to angular size of 5° x 5° and 
z direction scales between 0-1 in dimensionless units. The ampli- 
tude, A, of the GDSE is normalized i n the wa y described above to 
match 7^(325 MHz) = 20 K±2% dde Bruvn e t. al 1998), while B 
is normalized at 100 MHz: (5 = -2.55 ± 0.1 l lShaver et al]| 19991) . 

Fig. |2] shows a simulated map of the Galactic diffuse syn- 
chrotron emission according to the procedure described above, at 
a frequency of 120 MHz with an angular size of 5° x 5° on a 512^ 
grid. The mean brightness temperature of the map is Tt — 253K 
with a = 1.3K. 

In contrast to Fig.|2l which shows the angular variations of the 
GDSE at one frequency. Fig. [3] shows the amplitude variations of 
GDSE as a function of frequency for a number of lines of sight. 
Each line of sight has a slightly different power law index along 
the frequency direction as a result of the spatial variations in the 
temperature spectral index. Furthermore, the brightness tempera- 
ture variation for one line of sight is not a single power law but 
superposition of many power laws, due to the spectral index varia- 
tions both spatially and in the frequency direction. Note that Tt is 
still a very smooth function of frequency. 

3.2 Emission from SNRs 

Supernova remnants are composed of expanding shells that have 
strong magnetic fields which are able to produce cosmic rays. As 



Table 1. Angular size, flux density at 150 MHz (S150 MHz), spectral in- 
dex (q) and position on the map of the two s imulated supe rnova remnants. 
Values are calculated according to the data in lGreeJ )200d) catalogue. 
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Sl50 MHz 
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position on the map 




[arcmin] 


[Jy] 




[arcmin, arcmin] 


SNRi 


14 X 11 


7.91 


-0.65 


(254,53) 


SNRn 


5x6 


14.30 


-0.4 


(102,212) 



the particles escape the expanding shell, their energy decreases due 
to synchrotron cooling and we detect them at radio frequencies. The 
majority of the Galactic SNRs are within the Galactic plane but 
their distribution exp onentially decreases with distance from the 
Galac tic plane, z, (e.g. ICaswell & Lerchell979l :lx u. Zhang. &Had 
that is, ~ e ^. Moreover, due to the interaction of 
SNRs with the interstellar medium their radio surface brightness 
E decreases with an increase of their diameter D and with an 
increase of their height z, (e.g. ICaswell & Lerchjl 19791) , namely 
E ~ D-^e"". 

Our goal is to calculate the expected number of known SNRs 
within a LOFAR-EoR observational window at high Galactic lati- 
tudes, using the known number of observed radio SNRs from the 
iGreerJbOOel) catalogue and assuming that their distribution follows 
A'^ ~ e~^. On average, we obtain between one and two known 
SNRs in each 5° x 5° observational window. Given the extended 
nature of the SNRs we include two of them in each window in order 
to examine the influence of bright extended sources on the calibra- 
tion process and foreground removal. 

The simulated SNRs assume a power law spectrum: 

where S,^ is the flux density of a SNR at frequency So is its value 
at normalization frequency uq, and a is the spectral index. 

The simulated SNRs are placed randomly on the map and their 
angular size, flux den sity and spectral index are arbitrary chosen 
from the lGreenI ^20061) catalogue. The SNRs are added on the map 
as disks with uniform surface brightness. 

Properties of the two SNRs included in our foreground simu- 
lations are shown in Table[T] 

3.3 Diffuse free-free emission 

The diffuse thermal (free-free) emission contributes only ~ 1 
per cent of the total foregro unds within the fre quency range of 
the LOFAR-EoR experiment dShaveretalJI 19991) . It arises due to 
bremsstrahlung radiation in very diffuse ionized gas, with a total 
emission measure o f about 5 pc cm ~^ at high Galactic latitudes 
and Te = 8000 K ( lRevnolds|[l990h . This gas is optically thin at 
frequencies above a few MHz, so its spectrum is well determined 
and has a temperature spectral index of /3 = — 2.1. 

At high Galactic latitudes. Ha and free-free emission of 
the diffuse ionized gas are both proportional to the emission 
measure. Therefore, the Galactic Hq survey is generall y used 
as a tracer of the Galactic diffuse free-free emission jSmooll 
Il998l) . However, some groups also find significant correlation be- 
tween free-free emission an d dust emission ( Kogut et al. 199^; 
Ide Oliveira-Costa et alj[l997l) which can also be used as another 
independent tracer of the Galactic free- free emission. 

In our simulations we fo llowed iTegmark et alj ( l200(]h and 
ISantos. Coorav. & Knoxl ( l2005h who included the Galactic diffuse 
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free-free emission as a separate component of tlie Galactic fore- 
grounds witli an angular power spectrum Ci ~ and frequency 
Tb ~ v~^'^^ . Despite its small contribution to the foregrounds, the 
free-free emission is important for two reasons. Firstly, the ampli- 
tude of its angular fluctuations is much larger than that of the EoR 
signal. Secondly, and more importantly, its spectral index along 
the frequency direction is quite different from the other foreground 
components and could be important in testing the algorithms for 
the EoR signal extraction. 

To obtain the Galactic free-free emission maps we followed 
the same procedure as for the Galactic synchrotron emission with 
the additional simplification of fixing the power law index /3 to 
—2.15 across the map. is normalized according to the relation 
between Ha and free-free emission (see review bv lSmoot|[T99g) 
whereby the intensity of Ha emssion, la, is: 



7c = 0.36 R 



EM 



pc cm-6 J V 10* K 



(8) 



where EM is total emission measure and temperature. For < 
2.6 X 10'' K the value of 7 is 0.9. Combining Eq.[8]with the free- 
free equations in ISmootl l ll998l) .one finds a relation between and 



brightness temperature of free-free emission, Tj 



■//. 



r/^ (30 GHz) 



(9) 



Using Eq. [SjSEl together with EM = 5 pc cm~° and Te = 
8000 K, for high Galactic latitudes, one gets T/-''(30 GHz) = 
15.4 fiK. Assuming a frequency power law spectrum for Tj[^ with 
index -2.15, one obtains T/^' = 2.2 K at 120 MHz. 

Fig. |4] shows a simulated map of Galactic diffuse free-free 
emission at 120 MHz. The angular size of the map is 5° x 5° on 
512^^ grid, with the mean brightness temperature of Tt — 2.2K 
and a = 0.05K. 



4 EXTRAGALACTIC FOREGROUNDS 
4.1 Radio 



At the frequency range of the LOFAR-EoR experiment, bright 
radio sources are dominated by radio-loud galaxies, quasars and 
BL Lac objects (an AGN class of objects). However, at sub- 
mjy flux densities the contribution of late-type (star forming) 
galaxies, whose radio synchrotron emission originates from super- 
novae rather than AGN, becomes significant tPrandoni et al. 2001. : 
iMagliocchetti et al.ll2002l : ISadler et alj2002h 

The bright extragalactic radio sources are normally divided 
into two classes based on the relative physical position of their 
high and low surface brightness area within the lobes. These two 
classes are called FRI and FRII radio sources jFanaroff & RilevI 
ll974IJackson & Wall|l99"9h. O ur sim ulations of radio galaxies are 
based on the tables by Jackson ('2005') of extragalactic radio source 
counts at 151 MHz. Jackson (2005) has used ACDM based mod- 
els to calculate the evolution of the radio luminosity function of 
these sources, from which was predicted the source distributions 
and their number densities. 

In obtaining these tables. 'Jack soni ( [20051) assumed that the ra- 
dio sky consists of three population types of radio sources: FRI, 
FRII and star forming galaxies. Moreover, it is assumed that the 
local radio luminosity function of s tar forming galaxi es can be de- 
termined from the 2dFGRS-NVSS dSadler et al.i2002h galaxy sam- 
ple at 1.4 GHz. The parameterized number density and luminosity 
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Figure 4. Simulated map of total intensity emission of Galactic diffuse free- 
free emission with angular spectral index a = — 3 and frequency spectral 
index = —2.15. The angular size of the map is 5° X 5°, with ~ 0.6' 
resolution. The colour bar represents the brightness temperature Ti, of the 
Galactic diffuse free— free emission in Kelvin at 120 MHz. 



evolution of star forming galaxies is adopted from 'Haars ma et al] 
(2000). For t he local radio lu minosity function of FRI and FRn ra- 
dio galaxies, IjacksonI ( l2005l) assumed that it can be determined by 
exponential fitting of luminosity-dep endent density evolutio n to the 
observed source counts at 151 MHz d Jackson & Walll 19991). G iven 
the three evolving radio luminosity functions, 'jacksor/ (^2005^ sim- 
ulated the sky at different frequencies by randomly positioning 
and orientating each source on the sky. The intrinsic size is also 
selected randomly, assuming redshift independence. The FRI and 
FRII sources were modelled as double-lobe structures, and the star 
forming galaxies as circular discs. 

In our simula tions of radio ga laxies we adopt the three types of 
radio sources from |jacksonl ( l2005h and use the predicted source sur- 
face densities per deg'^ for 10, 5, 2, 1 and 0.1 mjy flux density limit 
in order to obtain the number of sources with certain flux density 
per deg^. However, and in contrast to the simulations by Jackson 
where each source is randomly positioned, we introduce a angular 
clusteri ng of the sources. The clustering i s motivated by the re- 
sults of iDi Matteo. Ciardi, & Miniatil ( |2004|) in which they showed 
that the contribution of the angular clustering of extragalactic ra- 
dio sources to the angular fluctuations of the foregrounds, at scales 
> 1 arcmin, is dominated by bright sources. Hence, in order to 
detect angular fluctuations in the cosmological 21-cm emission, ef- 
ficient source removal S > O.lmJy should be carried out. 

For angular clustering of the radio galaxies we used the partic- 
ul arly elegant procedure of Rayleigh-Levy random walk proposed 
bv lMandelbrol ( ll975Lll977l) . Starting from any arbitrary position, 
one places the next galaxy in a randomly chosen direction at angu- 
lar distance 0, drawn from the distribution: 
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Table 2. The flux density distribution of tlie simulated radio galaxies on the 
5° X 5° map (Fig. [5). The corresponding frequency is 150 MHz. 

number of sources with flux density limit 
10 mJy 5 mJy 2 mJy 1 mJy 100 /xJy 



>kFR1 oFR2 aSF 



FRi 20 
FRil 25 
SF 4 



55 122 
39 48 
38 210 



177 1001 
54 89 
419 18379 



{e/OoT for e-^Oo 

1 for 6* < 00, 7 > 0. 



(10) 



This is repeated many times until the correlation function of the 
distribution converges to the one desired. However, to compare the 
introduced correlation with obervational results and set the right 
values of 7 and Oq, the two point correlation function needs to be 
calculated. 

The two point correlation function, w{6), of the radio galaxy 
population is defined as the excess probability, over that expected 
for a Poisson distribution, o f finding a gala xy at an angular distance 
6 from a given galaxy (e.g. iPeebledl 19801) : 



&p = n[i + w{e)]m, 



(11) 



where 5P is probability, n the mean surface density and 50, a sur- 
face area element. Given the very large number of galaxies that can 
be simulated, we adopted the simplest form for estimating wiO), 
defined as: 



No{9) 



(12) 



where TVd is the number of pairs of galaxies with separation 6 in the 
correlated sample of galaxies and A'^r is the number of pairs with 
the same separation 6 but in a randomly distributed uncorrelated 
sample of galaxies. The total number of galaxies of the two samples 
must be the same. 

Recent results on the angular clustering of radio sources in the 
NRAO VLA Sky Survey (NVSS) ai id Faint Images of th e Radio 
Sky at Twenty centimetres (FIRST) jOverzier et alj[2003h showed 
that the two point correlation function is best fitted by a double 
power law w{e) = Z36l^"™ + AO'^"'-^ with slopes of 7b = 4.4, 
7A = 1.8 and amplitudes B = (1.5 ± 0.2) x 10"^ A = (1.0 ± 
0.2) X 10"^. However, in this study we adopt the simpler single 
power law correlation function with only two parameters, 7 = 1.8 
and A = 0.002. The reason for doing this is that steeper power 
law component (73) of the observed two point correlation function 
is dominant on angular scales smaller than those resolved by the 
LOFAR-EoR experiment. Note that 7 in the two point correlation 
function is the same one as in Eq.[TT] while A = (l/^o)"'- 

Fig- ID shows a simulated map of radio galaxies with angular 
power law distribution. All sources are point like as the angular res- 
olution of the LOFAR-EoR project will not be sufficient to resolve 
most of them. The angular extent of the map is 5° x 5°. On the 
map there are in total 20690 radio galaxies: 92 per cent SF radio 
galaxies, 7 per cent FRI and 1 per cent FRII radio galaxies. The 
flux density distribution at 150 MHz of the simulated radio galax- 
ies are shown in the Table[2l while the two point correlation fuction 
is w{e) = 0.0029-°-^ {A = 0.002, 7 = 1.8). The simulated ra- 
dio galaxies as sume a power l aw spectrum (see Eq.|7]( with spectral 
index -0.7 (e.g. |jacksonll2005h . 
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Figure 5. Simulated 2D angular clustering of radio galaxies (92 per cent 
- star forming (SF), 7 per cent - FRi and 1 per cent - FRii radio galaxies) 
with angular correlation function of the form w{9) = Q.0Q29~^'^ (A 

■ 7 = 1.8j. The flux density distribution of the galaxies is shown in 
the Tabled The angular size of the map is 5° X 5°. 



4.2 Radio clusters 

Galaxy clusters as radio sources are classified into cluster radio 
haloes and cluster radio relics. The former are morphologically reg- 
ular diffuse sources, typically centred inside the cluster and mostly 
unpolarized, whereas the radio relics are typically irregular, located 
at the periphery of the cluster and consist mostly of polarized ra- 
dio diffuse sources. Both types of clust er radio sourc e have steep 
frequency spectra with /9 ~ — 3 (see e.g. lFerettil2002l for review). 

The emission in radio haloes is due to synchrotron radiation 
by relativistic electrons with energies of ~ 10 GeV in mag- 
netic fields. The distribution of the radio haloes seems to follow 
closely the large sca le distribution of th e free-free driven X-ray 
emission of clusters ( iGovoni et ai]|200lh . This association is also 
supported by a strong correlation between th e radio halo luminos- 
ity and the host cluster X-ray luminosity (e.g.'EnBlin & RottgeringI 
I2OOZ) . However, not all clusters host radio haloes. Statistically, it 
is found that roughly 30-40 per cent of galaxy clusters with X-ray 
luminosity Lx 10*^ erg do host radio haloes. 

In our simulations of extragalactic foregrounds maps, as a 
starting point for simulating radio clusters, we used the ACDM 
deep wedge cluster catalogue of The Hubble Volume Projecfl The 
catalogue was obtained from an N-body simulation with one billion 
dark matter particles. This catalogue provides a list of clust ers up to 
redsh ift z ^ 4.37 with an gular coverage of 100 deg^ (Colbe rgetaU 



l200d : [Jenkins et al.l200ll : TEvrard et al j2003) . 

In order to translate the cluster mass into X-ray luminosity ix 



^ The cluster catalogue is part of simulations carried out by the Virgo 
Supercomputing Consortium using computers based at the Comput- 
ing Centre of the Max-Planck Society in Garching and at the Edin- 
burgh parallel Computing Centre. The data are publicly available at 
|http://www.mpa-garchiiig.mpg.de/galfonii/virgo/hubble| 
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and then into radio lumin osity Lr, we used the empiric al mass-X- 
ray luminosity relation o f iReiprich & Bohringeil j20o3) : 



Lx = axlO^^/igo^ erg s ^ 



Mr&b 



6x 



(13) 



and th e X-ray-radio luminosity relation of lEnfilin & Rottgerind 
( l2002h : 



ir,1.4 GHz = arlO^''/l50^ W Hz ^ 



(14) 



where ax = 0.449 6x = 1.9, flr = 2.78 and 6^ = 1.94 
tenfilin & Rottgerind [2OO2I) . It is important to note that Eq. [14] 
is derived for the radio frequency at 1.4 GHz. Note also that in 
Eq.[T3]the mass, Mr&b, is related to the cluster real mass, M, by 
AfR&B ~ MOU,^. 

Since Lr is derived for 1.4 GHz, we extrapolate the luminosity 
of each cluster to lower frequencies, relevant to the LOFAR-EoR 
experiment, according to: 



v\MHz\ 



(15) 



1400 

where a = -1.2 jKempner etar]|2004b . 

The angular size of the radio clusters was estimated from their 
physical radius and redshift. For the physical radius we used the 
virial radius r^jr calcu lated from the cluster mass according to 
( iBrvan & Normanll 19981) , 



(16) 



where Ac is the mean density and pcrit is the critical density at 
redshift z. 

In order to obtain the maps of radio clusters at the angu- 
lar and frequency ranges desired we first randomly choose 30 per 
cent of the catalogue's clusters (note that observations show that 
only 30-40 per cent of clusters have radio properties). The clus- 
ter masses are then used to estimate the radio luminosity of each 
cluster (Eq. |13|[T4l &|15b and its virial radius (Eq.|16t. Finally, the 
cluster is projected onto the simulated map according to its survey 
coordinates in the Hubble Volume simulation, its redshift and es- 
timated size. Note that the radio clsuters are added on the map as 
disks with uniform surface brightness. 

Fig. ID shows a 5° X 5° map of radio clusters simulated at 
120 MHz. The colour bar represents the brightness temperature of 
the clusters in logarithmic units. The size of each cluster has been 
scaled by a factor of 10 for visual clarity. 



5 POLARIZATION 

The need for understanding the polarization response stems from 
two main factors. One is the geometry of the LOFAR configura- 
tion, and the other is the cross-talk contamination between the two 
dipoles of a LOFAR antenna. In order to detect the EoR signal, 
which has at best a signal to noise ratio (SNR) of ~ 0.2 per beam 
on most angular scales - assuming 400 hours of observation un- 
der the configuration specified in Section [6]- we need to fully un- 
derstand the response of the LOFAR system in total intensity and 
polarization. As discussed before, this is vital for us in order to be 
able to span a dynamic range of 4—5 orders of magnitude. 

Since the LOFAR antennae are fixed to the ground, the sources 
are tracked only by beam-forming and not by steering the anten- 
nae as in traditional radio astronomy. This implies that, depend- 
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Figure 6. Simulated 120 MHz map of radio clusters of galaxies. The an- 
gular size of the map is 5° X 5°. Note that the colour bar represents the 
brightness temperature of the clusters in logarithmic units and that the size 
of each cluster has been scaled by a factor of 10 for visual clarity. We also 
assume uniform surface brightness of the clusters. 



ing on the position of the source on the sky, only a certain projec- 
tion of the two dipoles is visible. This projection changes as the 
source is tracked over time. Therefore, at most times the sources in 
the sky see different projections of each of these dipoles. Now, if 
the sources/foregrounds are polarized, we immediately see that the 
power output from the pair of dipoles (which is the sum of the two 
polarized components) will vary dramatically. This has to be taken 
in account almost exactly during the inversion and calibration pro- 
cesses in order to achieve the desired dynamic range. 

On the other hand, a leakage in the electronics can cause the 
power that is supposed to go through one of the dipoles to show up 
in the other (referred to as cross-talk). Although the cross-talk is 
small compared to the effect due to projection, we need to take it 
into account to eliminate any systematics. 

In this paper we show a first order simulation of the Galactic 
diffuse polarized emission and defer a more advanced discussion 
on the topic to a future paper. 

There are several polarization surveys of the Galactic syn- 
chrotron em ission between 327 MHz and 2.7 GHz (see review, 
iReicri I2OO6I) . The most recent one was done with the West- 
erbo rk telescope at 327 MHz, with arc minute angular resolu- 
tion jWieringa et al.l [19931 ; iHaverkom. Katgert, & de Bruvnilioool 
I2OO2I I2OO3I ) Its low frequency maps reveal a large amount of 
unusually shaped polarized small-scale structures, which have no 
counterpart in the total intensity. These structures are normally at- 
tributed to the coexistence of magnetic fields and thermal gas in the 
interstellar medium, which produce Faraday rotation at each line of 
sight. 
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The Faraday rotation depends on the observing frequency and 
rotation measure (RM) of the structure and it is defined along the 
line of sight. In order to measure Faraday rotation, observations at 
two or more frequencies are required. However, full understanding 
of the observed results could be quite difficult due to the possibil- 
ity of multiple Faraday rotation layers (screens) along t he lin e of 
sight and depolarization effects. ISrentiens & de Bruvnl ( 120051) in- 
troduced a new method (Faraday Rotation Measure Synthesis) that 
is able to cope with multiple screens and analyzes the contribution 
of each of these screens separately. 

The Galactic diffuse synchrotron emission is linearly polar- 
ized and its polarized intensity Ip can be expressed in terms of 
Stokes parameters U and Q: 



(17) 
(18) 



Ip = VU^ + 

. 1 u 
w = - arctan — 

where 9 is the polarization angle. 

In order to simulate the polarization maps of the Galactic 
synchrotron emission at low radio frequencies, we use a simple 
mod el of the Galactic syn chrotron polarization at high frequen- 
cies ( lGiardinoet"ai]|2002l) in combination with Faraday screens 
that are introduced to account for the effects of rotation and de- 
polarization at low frequencies. Note that the assumption about the 
correlation between the polarized and total intensity going into the 
Galactic synchrotron polarization model at high frequencies is not 
valid at low frequencies, since the observations mentioned above 
show polarized structures that have no counterpart in the total in- 
tensity. However, this assumption can be acceptable for a first es- 
timate. Note that the correlation assumption should not be used in 
the analysis of the redshifted 21 cm data as it is mostly invalid and 
can lead to wrong interpretations. In the future, we will improve on 
the model itself and use results from real polarization data obtained 
by the LFF^ 

The basic assumptions of the lGiardino et alj ( l2002h model of 
Galactic synchrotron polarization at high frequencies are: 

(i) The polarized component of Galactic synchrotron emission 
is proportional to the unpolarized intensity, which in terms of 
brightness temperature Tb is: 



Q = fn cos(2e) 

u = fn sm{2e) 



(19) 



where / is the fraction of polarized emission (or polarization de- 
gree) and 9 is the polarization angle; 

(ii) The fraction of po larized radiation / is related t o the tem- 
perature spectral index /3 jCortiglioni & Spoelstra|[l995l) : 



/ 



_ 3/3-3 
~ 3/3- 1 

(iii) The polarization angle 9 is given by: 



(20) 



9 =— arctan{x/ n,y/n) (21) 

where x, y are 2D random Gaussian fields with the mean zero and 
characterized by a power law spectrum, while n = y'x^ + y^. 
The power law spectral index is a — —1.7 and its value is 



° LFFE (Low Frequency Front End) are receivers at the Westerbork Syn- 
thesis Radio Telescope (WSRT) and cover the frequency range from 115 to 
170 MHz 
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Figure 7. Simulated 120 MHz map of polarized intensity (Ip) of diffuse 
Galactic synchrotron emission, and polarization angle (white lines). The 
angular size of the map is 5° X 5°, with ^ 0.6' resolution. 



driv en by observations, e.g. the Parkes 2.4 GHz polarimetric sur- 
vey jPuncan et all 199^ . 

The Faraday sc reens are modelled as 2D fi elds of rotation an- 
gles l\9 defined by dRvbicki &Lightmanlll986h : 



(22) 



where A is the wavelength of radiation and RM is the rotation mea- 
sure modelled as a 2D Gaussian random field (GRF) with a power 
law spectrum of spectral index a. Note that for the demonstrative 
purpose of this simulation we introduced only two Faraday screens, 
with mean zero and standard deviation 0.3, and arbitrarily set the 
value of a to -2. 

Therefore, in order to generate polarization maps of Galactic 
synchrotron emission at given frequencies, first we take the GDSE 
maps of total intensity (T) and temperature spectral index (/3) from 
Sec. 13. H and calculate the fraction of polarized radiation according 
to Eq. |20] Then, using Egs. I21l&l22lwe obtain polarization angle 
9 and Faraday rotation angle AS. Finally, we use Eq. [19] to get 
polarization maps Q and U. The angle in Eq.[T9]is the sum of A9 
over all Faraday screens and 9. 

Fig. |7] shows the simulated 120 MHz map of polarized inten- 
sity {Ip) of diffuse Galactic synchrotron emission. The polarization 
angles are shown as white lines. The Stokes Q map of simulated 
Galactic polarized emission is shown in Fig. [8] The related Stokes 
U map looks very similar to the Q map. 
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Figure 8. Simulated 120 MHz Stokes Q map of polarized diffuse Galactic 
synchrotron emission (DGSE). The angular size of the map is 5° X 5°, with 
~ 0.6' resolution. Simulated Stokes U map of polarized GDSE looks very 
similar to Q map. 



6 INSTRUMENTAL EFFECTS 

In this section we give a basic overview of the simulations of LO- 
FAR antenna response and show how the foreground maps are 
seen by LOFAR. More detailed discussion on the LOFAR response 
and data model fo r the LOFAR-EoR experiment will appear in 
iLabropoulos et all in prep.. For the LOFAR-EoR experiment we 
plan to use the LOFAR core, which will consist of approximately 
25 stations. However, in this paper we set the number of LOFAR 
core stations to 24. Each station is further split into two substa- 
tions which are separated by a few tens of metres (see Fig.|9}. Each 
substation consists of 24 tiles, with each tile having 4x4 crossed 
dipoles. For our goals we assume that each of the forty-eight sub- 
stations is a circular array with a radius of thirty-five metres. The 
stations are distributed in a randomized spiral layout and span a 
baseline coverage from 40 to 2000 metres. The total effective col- 
lecting area for the LOFAR-EoR experiment is ^ 0.07 km'^ at 
150 MHz. The instantaneous bandwidth of the LOFAR telescope is 
32 MHz and the aim for the LOFAR-EoR experiment is to observe 
in the frequency range between 115-180 MHz, which is twice the 
instantaneous bandwidth. To overcom e this, multiplexing in time 
has to be used (for more details see. lde Bruvn et. a]||2007f) . For the 
purpose of this paper we ignore this complication and assume 400 
hours of integration time for the hole frequency range. The LO- 
FAR specifications are not final yet and they might slightly change 
in the near future. A detailed description of the telescope together 
with the data model will be provided in forthcoming papers. 

In order to compute the true underlying visibilities, we make 
some simplifying assumptions. We assume that the narrow band- 
width condition holds and that the image plane effects have been 
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Figure 9. Position of 48 substations (24 stations) of the LOFAR core used 
for simulations of instrumental effects. Note that each substation consists of 
24 tiles, with each tile having 4x4 crossed dipoles. 



calibrated to a satisfactory level. This includes station complex gain 
calibration, a stable primary beam, and adequate compensation for 
the ionospheric effects, such that the ionospheric phase introduced 
during the propagation of electromagnetic waves in the ionosphere 
and the ionospheric Faraday rotation are corrected for. For an in- 
terferometer, the measured spatial correlation of the electric fields 
between two antermae is called 'visibi lity' and is approximately 
given bv lTavlor. Carilli & PerlevI h999t) : 



Vf{u,v)^ / A(Z,m)//(;,m)e'("'+""'^dWm 



where A is the primary beam, 1/ is the intensity map corresponding 
to frequency /, {u, v) are the coordinates, as seen from the source, 
of the tracks followed by an interferometer as the Earth rotates, and 
(/, m) are the direction cosines. 

We further treat each pixel of the map as a point source with 
the intensity corresponding to intensity of the pixel. Note that the 
equation above takes into account the sky curvature. The visibili- 
ties are sampled for all substation pairs and also at different pair 
positions, as the Earth rotates. 

We calculate the Fourier transform of the foreground model 
for each frequency in the above range. For every baseline and fre- 
quency, the uv tracks sample different scales of the Fourier trans- 
form of the sky at that frequency. Thus, the sampling function S 
becomes 



S{u,v) = ^ (5 (ll - Uk)S (v - Vk) , 



(23) 



where the summation is carried over all the pixels k. 

We compute those tracks for each interferometer pair for 4 
hours of synthesis with an averaging interval of 100 s and we then 
grid them on a regular grid in the uv plane. The maximum baseline 
assumed for the LOFAR core is 2 km and the station diameter is 
35 m, the number of independent elements in the uv plane is ~ 
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60'^. If the uv plane is oversampled by a factor of four, this yields 
256'^ pixels in the uv plane of « 60 rr? in size. After counting 
how many track points fall within each grid cell, we end up with 
a matrix representing the naturally weighted sampling function in 
the uv plane. By multiplying this sampling matrix with the Fourier 
transform of our model sky we get the visibilities on that grid with 
appropriate weights. This procedure is done for each baseline pair. 



Vf {u, v) = S-TI 



(24) 



where TI is the Fourier transform of the input image and S is the 
sampling function. 

The LOFAR visibility densities per resolution element at 
150 MHz, for the LOFAR-EoR experiment, are shown in Fig. 1101 
The total integration time is 400 hours (100 x 4 h night"^) with 
averaging time of 100 s and observing declinations 5 = 90° (left 
panel) and 5 — 52° (right panel). 

The inverse Fourier transform of the sampled visibilities is 
called the 'dirty' map. It is actually the sky map convolved with 
the Fourier transform of the sampling function, which is called the 
'dirty' beam or the 'PSF'. This is a simple-minded approach to 
estimating the sky brightness as it uses linear operations. The ap- 
proximation of the underlying brightness with the 'dirty' map is not 
always satisfactory, as side lobes from bright features will obscure 
fainter ones. In cases of low signal to noise, however - such as dur- 
ing the observation of the redshifted 21 -cm transition line of HI - 
one might choose not to proceed further than this first approxima- 
tion. To go beyond that we need extra information like the positivity 
of the intensity and compact support. The discussion of such issues 
is beyond the scope of this paper. This incomplete sampling of the 
uv plane also means that we do not measure the complete power at 
all scales, due to the holes in the uv coverage and its finite extent. 

An example of a 'dirty' map of the diffuse components in 
the foregrounds is shown in the Fig. [TT] together with the 'origi- 
nal' simulated foreground map with no interferometric effects and 
noise. The corresponding total integration time is 400 hours, with 
an averaging time of 100 s at 150 MHz. Note that the 'dirty' maps 
are generated without the inclusion of noise. 

The ultimate sensitivity of a receiving system is determined 
principally by the system noise. The discussion of the noise prop- 
erties of a complex receiving system like LOFAR can be lengthy, so 
we concentrate for our purposes on some basic principles. The the- 
oretical rms noise level in terms of flux density on the final image 
is given by 

_ 1 SEFD 

Vs X {N -1) X AV X tint 

where 77s is the system efficiency that accounts for electronic, dig- 
ital losses, N is the number of substations, Au is the frequency 
bandwidth and tint is the total integration time. SEFD is the System 
Equivalent Flux Density in Jy. The system noise we assume has 
two contributions. The first comes from the sky and is frequency 
dependent (~ i^^^''") and the second from receivers. 

For the LOFAR core the SEFD will be around 1000 Jy a t 
150 MHz, depending on the final design |de Bruvn et. a il l2007h . 
This means that we can reach a sensitivity of 520 mK at 150 MHz 
with 1 MHz bandwidth in one night of observations. In order 
to calculate the SEFD we use the following system temperature 
(Tsys) scaling relation as function of frequency (1/): Tays = 140 + 



^ This is the closest power of two to match the number of sampled ele- 
ments. By doing this one can benefit from the speed of the FFT. 



60(z^/300 MHz)~^'^^. Accumulating data from a hundred nights 
of observations brings the sensitivity down to 52 mK. We further 
assume that the distribution of noise over the map at one frequency 
is Gaussian. The noise contribution to each pixel in a map is drawn 
independently from a Gaussian distribution. The EoR signal is ex- 
tracted from two different scenarios. The first scenario involves the 
extraction of the signal from the 'original' maps - simulated maps 
that are not convolved with the dirty beam - after adding the noise 
directly to the 'original' maps. In the other scenario, the EoR sig- 
nal is extracted from 'dirty' maps to which we do not add noise but 
convolve the 'original' maps of the EoR signal plus the foregrounds 
with a simplified dirty beam. 

As the uv coverage scales linearly with frequency, one has to 
be careful in using the 'dirty' maps for extraction. This is because 
a pixel sampled at a given frequency need not be sampled at a later 
frequency. Since the analysis performed in this paper involves data 
across the frequency domain, we need a good uv coverage. If the 
uv sampling functions are not scaled accordingly, we will introduce 
additional difficulties arising from the mixing of spatial scales. To 
overcome spatial scale mixing, one of the strategies in the data anal- 
ysis is to use only the uv points that are present at all frequencies. 
In other words one can construct the uv plane mask that only con- 
tains the uv points that are sampled at every frequency. This step of 
course results in substantial data loss. 

The uv coverage for the LOFAR-EoR experiment changes in 
scale by ~ 40% between the frequency range of observation (115- 
180 MHz). By choosing only those uv points that are present at all 
frequencies, ~ 5% of the total data is lost in the frequency range 
specified above. Since with decreasing bandwidth of observation 
the amount of data lost decreases, one of the strategies could be to 
observe in windows of smaller bandwidth. However the observa- 
tional strategy of the LOFAR-EoR experiment is not yet finalized 
and will be discussed in detail in upcoming papers of the project. 
A detail discussion on the scaling of the uv coverage with fre- 
quency and its influence on the number of di scarded baselines a nd 
the amount of data loss will be discussed in iLabropoulos et all in 
prep.). 

In the following section we will show our ability to statisti- 
cally detect the EoR signal from the 'original' maps that include 
realistic level for the noise and from 'dirty' maps that do not in- 
clude the noise but are sampled with the uv mask that contains only 
uv points present at each frequency. In both cases the statistical de- 
tection of the signal is done on total intensity maps only. Moreover, 
perfect calibration is assumed and any other systematics that might 
influence the data are ignored. Those issues will be addressed in a 
follow-up paper. 



7 DETECTION OF THE EOR SIGNAL FROM THE FG 

This section presents the results on the statistical detection of the 
EoR signal from 'original' and 'dirty' LOFAR-EoR data maps that 
include the cosmological 21cm signal, diffuse components of the 
foregrounds and realistic noise. By 'original' maps we mean maps 
before inversion or in other words maps with no calibration errors 
or interferometric effects. 'Dirty' maps include only simplified uv 
coverage effect as an interferometric effect, but no calibration er- 
rors. 

By using only diffuse components of the foregrounds (Galac- 
tic diffuse synchrotron and free-free emission and integrated emis- 
sion from unresolved extragalactic sources) we assume that all re- 
solved discrete and extended sources have been successfully re- 
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Figure 10. The expected LOFAR visibility densities per resolution element at 150 MHz/or 400 h of total integration time f 100 X 4 h night ^ j with averaging 
time of 100 s and for observing declinations & = 90° (left panel) and 5 = 52° ( right panel). 
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Figure 11. Total intensity map of the simulated diffuse components of the foregrounds ('original' map with no interferometric effects and noise) and its 
corresponding 'dirty' map after 400 h of total integration time with averaging of 100 s at 150 MHz. 



moved from the observed maps, without any subtraction residuals. 
Also note that our analysis is done on total intensity maps only. The 
polarized case will be considered in the follow-up paper. 

The foreground and noise maps are simulated in the frequency 
range between 115 MHz and 178.5 MHz in steps of 0.5 MHz. The 
original maps simulated for a 5° x 5° field on a 512^ grid are re- 



binned to a 128^ grid, so that each pixel corresponds to 2.3' which 
is the resolution attained by the core of the LOFAR telescope. 

The EoR maps are simulated between the frequencies of 
115 MHz and 178.5 MHz in steps of 0.5 MHz, corresponding 
to redshifts between 1 1.5 and 6.5. 

The mean of the EoR signal, foreground and noise maps at 
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Figure 12. One line of sight (one pixel along frequency) of the 'original' 
LOFAR-EoR data maps (upper solid black line), smooth component of the 
foregrounds (dotted black line), fitted foregrounds (dashed green line) and 
residuals (lower solid black line) after taking out the foregrounds. Note that 
the residuals are the sum of the EoR signal and the noise. 



each frequency are set to zero since LOFAR is an interferometric 
instrument and measures only fluctuations around the mean. The 
typical variations, a, over the map at 150 MHz are ~ 5 mK for the 
EoR signal, ~ 2 K for the foregrounds and ~ 52 mK for noise. 
Hereafter, these values are considered fiducial values for the EoR 
signal, foregrounds and noise. 

The analysis on the LOFAR-EoR data maps can be done in 
two ways: firstly along the frequency direction where the fore- 
grounds are assumed to be smooth in contrast to the EoR sig- 
nal; and secondly in the spatial domain where the EoR signal and 
some components of the foregrounds are spatially correlated, but 
the noise is not. In this paper we will demonstrate statistical detec- 
tion of the signal by analysis along the frequency direction, taking 
lines of sight (map pixels) one by one. 

Fig.[T2]shows one line of sight (one pixel along frequency) of 
the 'original' LOFAR-EoR data cube (upper black solid line) with- 
out interferometric effects. The first step in the extraction of the 
EoR signal is to take out the smooth foreground component (dot- 
ted black line). It is important to note, however, that the smooth 
component of the foregrounds is not a simple power law but a 
superposition of three power laws (Galactic synchrotron and free- 
free emission and integrated emission from unresolved extragalac- 
tic sources) including the fact that one of the power law indices /3 
(Galactic synchrotron emission) changes slightly with frequency. 

The simplest method for foreground removal is a polynomial 
fitting in logarithmic scale (log(rf,) — log(i^)). The dashed green 
line on Fig. [12] represents the foregrounds fitted with a 3'^'^ order 
polynomial in the logarithmic scale. 

Fig.|13|shows a comparison between the detected and original 
EoR signal for randomly chosen lines of sight in the case of the 
fiducial foreground level and without noise for 'original' maps. As 
one can see, there is an almost exact agreement: this confirms that 
our approach when applied to noiseless data does not introduce any 
systematic biases. 
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Figure 13. Detection of the EoR signal from the simulated foreground maps 
('original' maps), without interferometric effects and noise: for a single 
random line of sight ( top panel) and as a standard deviation over all lines of 
.sight (bottom panel). The solid green line represents the original simulated 
EoR signal, and the dashed black line the extracted EoR signal. 



After taking out the fitted foregrounds from the 'original' data 
maps, the residuals should contain only the noise and the EoR sig- 
nal (lower solid black line on Fig. I12t. However, the assumption 
here is that we have fitted well enough such that the residuals be- 
tween the fitted and the 'real' foregrounds are smaller than the EoR 
signal. Otherwise the EoR signal could be fitted out if we are over- 
fitting, or be dominated by the foreground fitting residuals if we are 
under-fitting the foregrounds. 

The extraction of the EoR signal from the residuals along one 
line of sight is an impossible task, since the level of the noise is 
order of magnitudes larger than EoR signal and its value is un- 
known for a certain pixel. However, general statistical properties of 
the noise (standard deviation as a function of frequency) might be 
determined from the experiment and be used to statistically detect 
the EoR signal. By statistical detection we mean determination of 
the standard deviation of the EoR signal over the entire map as an 
excess variance over the variance of the noise. The general statis- 
tical properties of the noise might be determined in two ways. The 
first method is based on the difference between measured fluxes of 
a discrete point source, with well know properties, at two consecu- 
tive frequency channels. The second one is based on the difference 
between the measured flux in total and polarized intensity at the 
same frequency. However, the accuracy of the both methods need 
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Figure 14. Statistical detection of the EoR signal from the 'original' 
LOFAR-EoR data maps that include diffuse components of the foregrounds 
and realistic noise fcrnoise(150 MHz) = 52 mKJ but without interfero- 
metric effects. The dashed-dotted black line represents the standard devi- 
ation (cr) of the noise as a function of frequency, the dotted green line the 
a of the residuals after taking out the smooth foreground component, and 
the solid red line the a of the original EoR signal. The grey shaded zone 
shows the la detection, whereas the dashed white line shows the mean of 
the detection. Note that the y-axis is in logarithmic scale. 



to be tested for the LOFAR-EoR experiment and we leave further 
discussion on this topic for a forthcoming paper. 

Fig. [14] shows the standard deviation of residuals as a func- 
tion of frequency (dotted green line), after taking out the smooth 
component of the foregrounds, by polynomial fitting in logarith- 
mic scale to each line of sight of the 'original' maps. The most 
satisfactory result we get with a order polynomial. The dashed- 
dotted black line represents the standard deviation of the noise. By 
subtracting (in quadrature) the anoisc from (Trcsiduais, we get the 
excess variance (cteor) of the EoR signal. However, in order to 
determine the error on the detection of the EoR signal, we con- 
ducted a Monte-Carlo simulation of the extraction of the signal. 
We made 1000 independent noise cube realisations and applied the 
signal extraction algorithm to each. The results of the simulation 
are shown in Fig. [14] The grey shaded zone shows the 2a detec- 
tion, whereas the dashed white line shows the mean of the detec- 
tion. As one can see the mean of the detected EoR signal is in good 
agreement with the original (solid red line) up to 165 MHz. The dis- 
agreement for higher frequencies is due to over-fitting and low EoR 
signal level. Remember that for most of the sightlines our simulated 
Universe has already been ionized at this frequency, corresponding 
to 2 ~ 7.5 (see section|2j. 

In order to see the influence of the foreground and noise level 
on the EoR extraction and detection scheme, we repeated the same 
analysis on 'original' maps of the four different models. The first 
model has a foreground level two times bigger than fiducial and 
the second four times; the third has the fiducial foreground level 
but smaller noise level by \/2; and the last one has a normal fore- 
ground level and no noise (see Table[3]l. Note that by fiducial fore- 
ground level and noise level we mean (7/g(150 MHz) ~ 2 K and 



Table 3. Five different sets of values for standard deviation of foregrounds 
{afg[K] at 150 MHz) and of noise (anoise[mK] at 150 MHz), used for 
testing the EoR extraction and detection scheme. Note that case (a) repre- 
sents the fiducial case. 





case (a) 


case (b) 


case (c) case (d) case (e) 




2 


4 


8 2 2 


o-„oise[mK] 


52 


52 


52 36 



o"noiBG(150 MHz) ~ 52 mK. The results are shown in Figs.ll3l& 

m 

Comparing Figs. [14] & [15] we see the higher foreground levels 
decrease the quality of the EoR detection. Lower quality in the EoR 
detection is due to over-fitting at higher frequencies. However, even 
for the four times bigger foreground level we are able to detect the 
EoR signal up to 150 MHz. 

Comparing Figs.[T3][T4]c&[T5] we can see that a lower noise 
level increases the quality of the EoR detection, as expected. Better 
precision in the EoR detection with lower noise level also confirms 
that our foreground removal procedure works well. 

Finally, in Fig. [TT] we show the statistical detection of the 
EoR signal from the 'dirty' foregrounds -I- EoR signal maps without 
any noise. Note that the 'dirty' maps are produced with a sampling 
function (uv mask) that contains only uv points present at each fre- 
quency, in order to overcome additional difficulties from the mixing 
of angular scales in the frequency direction introduced by the linear 
frequency variation of the uv coverage and incomplete sampling in 
the frequency direction. 

The smooth component of the foregrounds is removed by 
polynomial fitting to each line of sight. The most satisfactory result 
we get with a 6* order polynomial. A different order of polynomial 
from the case of the 'original' maps is required due to the angular 
scale mixing over each map introduced by convolution of the map 
with a 'dirty' beam. 

Fig. [I6] shows a comparison along the frequency direction of 
the same pixel from the 'original' (solid line) and 'dirty' (dashed 
line) foregrounds -I- signal maps. Note that the foregrounds are still 
smooth along the frequency direction of the 'dirty' maps, but the 
shape of the function is slightly different. The difference is due to 
incomplete uv coverage sampling and its finite extent, determined 
by the shortest and longest baselines. 

The dashed green line in Fig. [TT] shows the standard devia- 
tion, as a function of frequency, of the extracted EoR signal from 
'dirty' foregrounds -I- EoR signal maps. The dashed-dotted black 
line shows the standard deviation of the 'original' EoR signal maps, 
while the solid black line shows the standard deviation of the 
'dirty' EoR signal maps. The agreement between the standard devi- 
ations of the extracted and 'dirty' EoR signals is satisfactory, while 
their slightly lower levels than the 'original' signal are due to the 
smoothing effect of the instrumental response. 



8 DISCUSSION AND OUTLOOK 

This paper presents foreground simulations tailored for the 
LOFAR-EoR experiment that is set to study the redshifted 21 -cm 
hyperfine line of neutral hydrogen from the Epoch of Reionization. 
The foreground simulations include Galactic diffuse synchrotron 
and free-free emission, synchrotron emission from Galactic super- 
nova remnants and extragalactic emission from radio galaxies and 
clusters. For each of the foreground components, we generate the 
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Figure 15. Statistical detection of the EoR signal from the 'original' LOFAR-EoR data maps with foreground level two {left panel) and four {middle panel) 
times bigger than the fiducial foreground level, and with noise level smaller by \/2 f right panel) than the fiducial noise level, hut without interferometric effects. 
Colours and line coding are the same as in Fig. \l4\ Note that the y-axis is in logarithmic scale. 




Figure 16. One line of sight {one pixel along frequency) of the 'dirty' fore- 
ground {EG) + cosmological 21-cm signal (CS) maps (dashed line) in com- 
parison with the same pixel along the frequency of the 'original' EG+CS 
maps (solid line). The difference between these two lines is due to incom- 
plete uv coverage and its finite extent. 



Figure 17. Detection of the EoR signal from the simulated simplified 'dirty ' 
foreground maps, without noise, as a standard deviation a over all lines 
of sight (dashed green line). The dashed-dotted black line represents the a 
over all lines of sight of 'original' EoR signal, while the solid black line the 
a over all lines of sight of 'dirty' EoR signal maps. 



5° X 5° field in the frequency range approximately between 115 
and 180 MHz pertaining to the LOFAR-EoR. 

Since the diffuse Galactic synchrotron emission is the domi- 
nant component (~ 70 per cent) we include all its observed charac- 
teristics: spatial and frequency variations of brightness temperature 
and its spectral index, and brightness temperature variations along 
the line of sight. Discrete sources of Galactic synchrotron emission 
are included as observed emission from supernovae remnants. 

Despite the minor contribution of the Galactic free-free emis- 
sion (~ 1 per cent), it is included in our simulations of the fore- 
grounds as an individual component. It has a different temperature 
spectral index from Galactic synchrotron emission. 

Integrated emission from extragalactic sources is decomposed 
into two components: emission from radio galaxies and from ra- 



dio clusters. Simulations of radio galaxie s are based on the source 
count functions at low radio frequency bv ljacksodJiooi) . for three 
different types of radio galaxies, namely FRI, FRII and star forming 
galaxies. Correlations obtained by radio galaxy surveys are used for 
their angular distribution. Simulations of radio clusters are based on 
a cluster catalogue from the Virgo consortium and observed mass- 
X-ray luminosity and X-ray-radio luminosity relations. 

Under the assumption of perfect calibration, LOFAR-EoR 
data maps that include the simulated cosmological 21cm sig- 
nal ((7eor(150 MHz) ~ 5 mK), diffuse components of the 
foregrounds (ctfgs (150 MHz) ~ 2 K) and realistic noise 
(o"noisc(150 MHz) ~ 52 mK) are produced. We refer to this set of 
parameters as our fiducial model. For noise we assume it has two 
components, the sky noise and receiver noise. The former varies 
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with frequency as v^'^'^^ whereas the latter is roughly frequency 
independent. 

The extraction of the EoR signal is performed along the fre- 
quency direction, taking lines of sight (map pixels) one by one. 
The first step in the EoR extraction is removal of the smooth fore- 
grounds component for each line of sight (see Fig.ll2t. In our anal- 
ysis we fit a 3"* order polynomial in the logarithmic scale. How- 
ever, one should be careful in choosing the order of the polynomial 
to perform the fitting. If the order of the polynomial is too small, 
the foregrounds will be under-fitted and the EoR signal could be 
dominated and corrupted by the fitting residuals, while if the order 
of the polynomial is too big, the EoR signal could be fitted out. 

After foreground removal, the residuals are dominated by in- 
strumental noise. Since the noise is unknown for each line of sight 
and is an order of magnitude larger than the EoR signal, it is an 
impossible task to directly extract the EoR signal for each line of 
sight. However, assuming that the statistical properties of the noise 
(cTnoisc) will be known, we can use it to statistically detect the EoR 
signal. The statistical detection of the EoR signal is the measure of 
the excess variance over the entire map, (TeoRi that should be ob- 
tained by subtracting the variance of the noise, (J^aisc from that of 
the residuals, cr^^csiduais- 

Fig-HHshows the results of a successful statistical detection of 
the EoR signal in the fiducial model of the foregrounds and noise. 
The detected standard deviation of the signal is in a good agree- 
ment with original signal up to 165 MHz. The disagreement for 
higher frequencies is due to over-fitting caused by the very weak 
cosmological signal at these frequencies. 

In order to see the influence of the foreground and noise level 
on the EoR extraction and detection, the same analysis was done 
for models with two and four times bigger foreground levels than 
in the fiducial model, for a model where the noise is smaller by \/2, 
and for a model without noise (see Tabel[3}. The results are shown 
in Figs.[T5]&[l3] 

In the case of higher foreground levels, the EoR signal de- 
tection is hampered by over-fitting. In the case of lower noise lev- 
els, however, the proposed EoR detection algorithm performs ex- 
tremely well. 

For the diffuse components of simulated foregrounds, a 'dirty' 
map with realistic but idealised instrumental response of LOFAR is 
produced (see Fig. II It. However, the signal extraction scheme we 
apply only to the 'dirty' maps that have been produced with a uni- 
form uv coverage as a function of frequency and have no noise. This 
is due to the additional difficulties introduced by mixing of angular 
scales in the frequency direction. Those issues will be discussed in 
a follow-up paper. 

In addition to the simulations of the total brightness temper- 
ature, the polarized Galactic synchrotron emission maps are pro- 
duced. Here, we follow a simple model that includes multiple Fara- 
day screens along the line of sight (see Figs.|7]&[lll. The motivation 
behind these simulations is that improper polarization calibration 
could severely contaminate the EoR signal, so future robust extrac- 
tion algorithms have to take this into account. 

Fig.[T8]shows the angular power spectra of the simulated EoR 
signal (dotted red line), simulated diffuse component of the fore- 
grounds (solid black line) and three levels of noise (blue lines) at 
150 MHz. The lines are drawn as the best fit to the corespond- 
ing points. The dashed blue line represents the level of the noise 
(o-noise(150 MHz) = 52 mK) after one year of the LOFAR-EoR 
experiment (400 h of total observing time) with a single beam. 
For this case of instrumental noise and inclusion of realistic dif- 
fuse foregrounds we have shown that we are able to statistically 
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Figure 18. Angular power spectra of the simulated EoR signal (dotted red 
line), simulated dominant component of the foregrounds (solid black line) 
and three levels of noise: the dashed blue line represents noise for a sin- 
gle beam after one year of integration, the dashed-dotted blue line for five 
beams after one year of integration and the dashed-dotted-dotted blue line 
for five beams and four years of integration. The lines are drawn as the best 
fit to the coresponding points. 

detect the EoR signal despite the small signal to noise ratio. How- 
ever, the current obse rving plan of the LOFAR-EoR experiment 
lie Bruvn etraill2007h is to observe with five independent beams, 
which reduces the (Tnoisc by a factor of ^/Z (dashed-dotted blue 
line). After four years of observations (4 x 400 h) with five beams 
the (TnoisG is reduced by a factor of a/20 (dashed-dotted-dotted blue 
line), which means that the signal to noise ratio is roughly 0.5. 

Finally we would like to emphasize that this paper is just 
a first step in testing and developing the analysis tools for the 
LOFAR-EoR experiment. Future papers will use the foregrounds 
simulations deve loped in this p aper together with simulations of 
the EoR signal {Thomas et al.L in prep.), instrumental response 
jLabropoulos et all in prep.), and other non-astrophysical effects 
(e.g. ionosphere, RFIs, ...) in order to test all aspects of the pipeline 
in the LOFAR-EoR experiment. 
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